! subroutine dgesv (integer N, integer NRHS, double precision, dimension(
  !     lda, * ) A, integer LDA, integer, dimension( * ) IPIV, double
   !    precision, dimension( ldb, * ) B, integer LDB, integer INFO)
    !    DGESV computes the solution to system of linear equations A * X = B
     !  for GE matrices

module interface_d
  interface
     subroutine read_matrix(A, B)
       implicit none
        integer, parameter :: dbl = 8
        real(dbl), allocatable, dimension(:,:), intent(inout) :: A
        real(dbl), allocatable, dimension(:,:), intent(inout) :: B
      end subroutine read_matrix
 
      subroutine print_matrix(A)
       integer, parameter :: dbl = 8
       real(dbl), dimension(:,:), intent(in) :: A
     end subroutine print_matrix
  end interface
end module interface_d


program dgesv_
  use interface_d
  implicit none
  integer, parameter :: dbl = 8, fnum = 4
  real(dbl), allocatable, dimension(:,:) :: A, B
  integer :: n, nrhs, lda, ldb, info
  integer, allocatable, dimension(:) :: ipiv
  call read_matrix(A, B)
  call print_matrix(A)
  call print_matrix(B)

  n = size(A,1)
  nrhs = size(B,2)
  ldb = n; lda = n
  allocate(ipiv(n))
  call dgesv(n, nrhs, A, lda, ipiv, B, ldb, info)
  if(info == 0) then
     write(*,*) 'ok'
  end if
  call print_matrix(B)
  deallocate(A)
  deallocate(B)
  deallocate(ipiv)
end program dgesv_


subroutine read_matrix(A, B)
  implicit none

  integer, parameter :: dbl = 8, fnum = 4
  integer :: dim, b_col
  real(dbl), allocatable, dimension(:,:), intent(inout) :: A
  real(dbl), allocatable, dimension(:,:), intent(inout) :: B
  integer :: k, status

  open(unit=fnum, file='matrix.txt', status='old', action='read', iostat=status)
  openif: if(status /= 0) then
     write(*,*) 'open file error'
     stop
  else
     read(fnum,*) dim, b_col
     write(*,*) dim, b_col
     allocate(A(dim,dim))
     allocate(B(dim,b_col))
     
     readA: do k = 1,dim
        read(fnum,*) A(k,:)
     end do readA
 
     readB: do k = 1, b_col
        read(fnum,*) B(:,k)
     end do readB
  end if openif

  close(fnum)
end subroutine read_matrix

! 按行打印
subroutine print_matrix(A)
  integer, parameter :: dbl = 8
  real(dbl), dimension(:,:), intent(in) :: A
  integer :: row, k
  row = size(A,1)
  do k = 1, row
     write(*,*) A(k,:)
  end do
  write(*,*)
end subroutine print_matrix
